// Copyright 2022 Jeroen van Nugteren

// Permission is hereby granted, free of charge, to any person obtaining a copy of this software
// and associated documentation files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use, copy, modify, merge, publish,
// distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included in all copies or
// substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
// OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
// CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

// include header
#include "distmesh2d.hh"

#include "rat/common/elements.hh"

// code specific to Rat
namespace rat{namespace dm{

	// constructor
	DistMesh2D::DistMesh2D(){
		set_h0(RAT_CONST(0.2));

		// default scaling fun (no scaling)
		fh_ = DFOnes::create();
	}

	// factory
	ShDistMesh2DPr DistMesh2D::create(){
		return std::make_shared<DistMesh2D>();
	}

	// meshgrid creates an x,y grid from input arrays
	// same as meshgrid in matlab
	void DistMesh2D::meshgrid(arma::Mat<fltp> &x, arma::Mat<fltp> &y, 
		const arma::Row<fltp> &xgv, const arma::Col<fltp> &ygv){

		// allocate coordinates
		x.set_size(ygv.n_elem,xgv.n_elem);
		y.set_size(ygv.n_elem,xgv.n_elem);

		// set coordinates to matrix
		x.each_row() = xgv;
		y.each_col() = ygv;
	}

	// circumcircle creation
	// translated from from Paul Bourke 
	// http://paulbourke.net/papers/triangulate/
	bool DistMesh2D::circumcircle(
		arma::Row<fltp>::fixed<2> &pc, fltp &r,
		const arma::Row<fltp>::fixed<2> &p, 
		const arma::Mat<fltp>::fixed<3,2> &ptri){

		//  get points
		const fltp x1 = ptri(0,0);
		const fltp y1 = ptri(0,1);
		const fltp x2 = ptri(1,0);
		const fltp y2 = ptri(1,1);
		const fltp x3 = ptri(2,0);
		const fltp y3 = ptri(2,1);

		// allocate circle centers
		fltp xc,yc;
		fltp eps = RAT_CONST(1e-7);

		// check for coincident points
		if(std::abs(y1 - y2) < eps && std::abs(y2 - y3) < eps)
			return(false);

		// not coincident
		if(std::abs(y2-y1) < eps){ 
			const fltp m2 = - (x3 - x2) / (y3 - y2);
			const fltp mx2 = (x2 + x3) / RAT_CONST(2.0);
			const fltp my2 = (y2 + y3) / RAT_CONST(2.0);
			xc = (x2 + x1) / RAT_CONST(2.0);
			yc = m2 * (xc - mx2) + my2;
		}else if(std::abs(y3 - y2) < eps){ 
			const fltp m1 = - (x2 - x1) / (y2 - y1);
			const fltp mx1 = (x1 + x2) / RAT_CONST(2.0);
			const fltp my1 = (y1 + y2) / RAT_CONST(2.0);
			xc = (x3 + x2) / RAT_CONST(2.0);
			yc = m1 * (xc - mx1) + my1;
		}else{
			const fltp m1 = - (x2 - x1) / (y2 - y1); 
			const fltp m2 = - (x3 - x2) / (y3 - y2); 
			const fltp mx1 = (x1 + x2) / RAT_CONST(2.0); 
			const fltp mx2 = (x2 + x3) / RAT_CONST(2.0);
			const fltp my1 = (y1 + y2) / RAT_CONST(2.0);
			const fltp my2 = (y2 + y3) / RAT_CONST(2.0);
			xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2); 
			yc = m1 * (xc - mx1) + my1; 
		}
		
		// set circle center point
		pc(0) = xc;	pc(1) = yc;

		// check 
		fltp dx = x2 - pc(0);
		fltp dy = y2 - pc(1);
		fltp rsqr = dx*dx+dy*dy;
		r = std::sqrt(rsqr); 
		dx = p(0) - xc;
		dy = p(1) - yc;
		fltp drsqr = dx*dx + dy*dy;
		// return((drsqr <= rsqr) ? true : false);
		return((drsqr - rsqr) <= eps ? true : false);
	}


	// delaunay triangulation
	// translated from from Paul Bourke 
	// http://paulbourke.net/papers/triangulate/
	arma::Mat<arma::uword> DistMesh2D::triangulation(const arma::Mat<fltp> &p){
		// check input
		assert(p.is_finite());
		assert(p.n_cols==2);

		// settings
		arma::uword emax = 200; // edge buffer

		// get number of points
		arma::uword num_points = p.n_rows;
		arma::uword num_tri_max = 4*num_points;

		// allocate output
		arma::Mat<arma::uword> tri(num_tri_max,3);
		arma::uword num_tri = 0;

		// allocate edges
		arma::Mat<arma::sword> edges(emax,2);
		arma::uword num_edges = 0;

		// Find the maximum and minimum vertex bounds.
		// This is to allow calculation of the bounding triangle
		arma::Row<fltp> pmax = arma::max(p,0);
		arma::Row<fltp> pmin = arma::min(p,0);

		// size of box
		arma::Row<fltp> dp = pmax-pmin;

		// find maximum dimension
		fltp dmax = arma::max(dp);

		// midpoint
		arma::Row<fltp> pmid = (pmax + pmin)/2;

		// Set up the supertriangle
		// his is a triangle which encompasses all the sample points.
		// The supertriangle coordinates are added to the end of the
		// vertex list. The supertriangle is the first triangle in
		// the triangle list.
		arma::Mat<fltp>::fixed<3,2> psuper;
		psuper(0,0) = pmid(0) - 20*dmax;
		psuper(0,1) = pmid(1) - dmax;
		psuper(1,0) = pmid(0);
		psuper(1,1) = pmid(1) + 20*dmax;
		psuper(2,0) = pmid(0) + 20*dmax;
		psuper(2,1) = pmid(1) - dmax;

		// sort indexes for p
		arma::Col<arma::uword> sort_index = arma::sort_index(p.col(0));

		// combine points
		arma::Mat<fltp> pext = arma::join_vert(p.rows(sort_index),psuper);

		// set first triangle
		tri(0,0) = num_points;
		tri(0,1) = num_points+1;
		tri(0,2) = num_points+2;
		num_tri = 1;

		// keep track of completed triangles
		arma::Col<arma::uword> is_complete(num_tri_max);
		is_complete(0) = false;

		// include each point one at a time into the existing mesh
		for(arma::uword i=0;i<num_points;i++){
			// get point from list
			arma::Row<fltp>::fixed<2> pnt = pext.row(i);
			num_edges = 0;

			// Set up the edge buffer.
			// If the point (xp,yp) lies inside the circumcircle then the
			// three edges of that triangle are added to the edge buffer
			// and that triangle is removed.
			// walk over triangles
			for(arma::uword j=0;j<num_tri;j++){
				// check if is already complete
				if(is_complete(j))continue;

				// get points for this triangle from list
				arma::Mat<fltp>::fixed<3,2> ptri = pext.rows(tri.row(j));

				// check if point is inside circle
				arma::Row<fltp>::fixed<2> pc = {0,0}; fltp rc = 0;
				bool inside = circumcircle(pc,rc,pnt,ptri);

				// check if complete
				if (pc(0) + rc < pnt(0)){
					is_complete(j)=true;
				}

				// if point is inside circle
				if(inside){
					// check if edge list full
					if(num_edges+3>=emax){emax += 100; edges.resize(emax,2);}
					
					// add new edges to list
					edges(num_edges+0,0) = tri(j,0);
					edges(num_edges+0,1) = tri(j,1);
					edges(num_edges+1,0) = tri(j,1);
					edges(num_edges+1,1) = tri(j,2);
					edges(num_edges+2,0) = tri(j,2);
					edges(num_edges+2,1) = tri(j,0);
					num_edges+=3;

					// copy last triangle and reduce list
					tri.row(j) = tri.row(num_tri-1);
					is_complete(j) = is_complete(num_tri-1);
					num_tri--; j--;
				}
			}


			// Tag multiple edges
			// Note: if all triangles are specified anticlockwise then all
			// interior edges are opposite pointing in direction.
			for(arma::uword j=0;j<num_edges-1;j++){
				for(arma::uword k = j+1;k<num_edges;k++){
					if((edges(j,0) == edges(k,1)) && (edges(j,1) == edges(k,0))){
						edges.row(j).fill(-1); edges.row(k).fill(-1);
					}

					// Shouldn't need the following, see note above
					else if((edges(j,0) == edges(k,0)) && (edges(j,1) == edges(k,1))){
						edges.row(j).fill(-1); edges.row(k).fill(-1);
					}
				}
			}

			// Form new triangles for the current point
			// Skipping over any tagged edges.
			// All edges are arranged in clockwise order.
			for(arma::uword j=0;j<num_edges;j++) {
				 // check if tagged
				 if(arma::any(edges.row(j)<0))continue;

				 // add edges to new triangle
				 tri(num_tri,0) = edges(j,0);
				 tri(num_tri,1) = edges(j,1);
				 tri(num_tri,2) = i;
				 is_complete(num_tri) = false;
				 num_tri++;
			}
		}


		// Remove triangles with supertriangle vertices
		// These are triangles which have a vertex number greater than nv
		for(arma::uword i=0;i<num_tri;i++){
			if(arma::any(tri.row(i)>=num_points)){
				tri.row(i) = tri.row(num_tri-1);
				num_tri--; i--;
			}
		}

		// change size
		tri.resize(num_tri,3);

		// get unique rows
		// tri = cmn::Extra::unique_rows(tri);

		// unsort tri
		arma::Col<arma::uword> original_index = arma::regspace<arma::Col<arma::uword> >(0,num_points-1);
		original_index = original_index.rows(sort_index);
		tri = arma::reshape(original_index.rows(arma::vectorise(tri)),num_tri,3);

		// return triangulation
		return tri;
	}

	// element graph
	void DistMesh2D::element_graph(
		arma::Col<arma::uword> &row, 
		arma::Col<arma::uword> &col, 
		const arma::Mat<arma::uword> &t){

		// number of elements
		const arma::uword num_tri = t.n_rows;

		// Interior bars duplicated
		arma::Mat<arma::uword> e(3*num_tri,2);
		arma::Col<arma::uword> id(3*num_tri);

		// create edges
		e.rows(0*num_tri,1*num_tri-1) = t.cols(arma::Row<arma::uword>{0,1});
		e.rows(1*num_tri,2*num_tri-1) = t.cols(arma::Row<arma::uword>{0,2});
		e.rows(2*num_tri,3*num_tri-1) = t.cols(arma::Row<arma::uword>{1,2});

		// keep track of element indices
		id.rows(0*num_tri,1*num_tri-1) = arma::regspace<arma::Col<arma::uword> >(0,num_tri-1);
		id.rows(1*num_tri,2*num_tri-1) = arma::regspace<arma::Col<arma::uword> >(0,num_tri-1);
		id.rows(2*num_tri,3*num_tri-1) = arma::regspace<arma::Col<arma::uword> >(0,num_tri-1);

		// sort edges
		arma::Mat<arma::uword> es = arma::sort(e,"ascend",1);

		// sort rows
		for(arma::uword i=0;i<es.n_cols;i++){
			const arma::Col<arma::uword> sort_index = arma::stable_sort_index(es.col(i));
			es = es.rows(sort_index); id = id.rows(sort_index);
		} 

		// allocate graph
		row.set_size(3*num_tri); col.set_size(3*num_tri);

		// walk over edges
		arma::uword cnt = 0;
		for(arma::uword i=1;i<es.n_rows;i++){
			// same edge
			if(arma::all(es.row(i)==es.row(i-1))){
				//row(cnt) = id(i); col(cnt) = id(i-1); cnt++;
				row(cnt) = id(i-1); col(cnt) = id(i); cnt++;
			}
		}

		// resize
		row.resize(cnt); 
		col.resize(cnt);
	}

	// calculate element quality
	arma::Col<fltp> DistMesh2D::element2quality(
		const arma::Mat<arma::uword> &tq, const arma::Mat<fltp> &p){
		// number of elements
		const arma::uword num_elements = tq.n_rows;
		const arma::uword num_sides = tq.n_cols;

		// angle of ideal polygon
		fltp poly_angle = (num_sides-2)*arma::Datum<fltp>::pi;
		fltp corner_angle = poly_angle/num_sides;

		// allocate
		arma::Col<fltp> Q(num_elements);

		// walk over elements
		for(arma::uword i=0;i<num_elements;i++){
			// allocate corners
			arma::Row<fltp> theta(num_sides);

			// calculate corners
			for(arma::uword j=0;j<num_sides;j++){
				// get edge vectors
				const arma::Row<fltp>::fixed<2> v1 = p.row(tq(i,(j+1)%num_sides)) - p.row(tq(i,(j+0)%num_sides));
				const arma::Row<fltp>::fixed<2> v2 = p.row(tq(i,(j+2)%num_sides)) - p.row(tq(i,(j+1)%num_sides));
				
				// dot product
				fltp costheta = -arma::sum(v1%v2)/(std::sqrt(arma::sum(v1%v1))*std::sqrt(arma::sum(v2%v2)));
				
				// fix bounds
				if(costheta<-1)costheta = -1;
				if(costheta>1)costheta = 1;

				// calcualte angle
				theta(j) = std::acos(costheta);
			}

			// concave elements
			if(std::abs(arma::sum(theta)-poly_angle)>=1e-4){
				Q(i) = -RAT_CONST(1.0);
			}else{
				Q(i) = arma::prod(RAT_CONST(1.0) - arma::abs((corner_angle - theta)/corner_angle));
			}
		}

		// return quality
		return Q;
	}


	// calculate element quality
	arma::Col<fltp> DistMesh2D::element2quality_alt(
		const arma::Mat<arma::uword> &tq, const arma::Mat<fltp> &p){
		// get two common points
		const arma::Mat<fltp> p1 = p.rows(tq.col(1));
		const arma::Mat<fltp> p3 = p.rows(tq.col(3));

		// calculate cross vector
		arma::Mat<fltp> v = p3-p1;

		// normalize
		v.each_col()/=arma::sqrt(arma::sum(v%v,1));

		// check direction
		const arma::Row<fltp>::fixed<2> t1 = {0,1};
		const arma::Row<fltp>::fixed<2> t2 = {1,0};

		// dot product
		const arma::Col<fltp> vt1 = arma::sum(v.each_row()%t1,1);
		const arma::Col<fltp> vt2 = arma::sum(v.each_row()%t2,1);

		// allocate
		const arma::Col<fltp> Q = RAT_CONST(1.0) - arma::max(arma::abs(vt1),arma::abs(vt2));

		// return quality
		return Q;
	}

	// merge triangles
	arma::Row<arma::uword>::fixed<4> DistMesh2D::merge2triangles(
		const arma::Row<arma::uword>::fixed<3> &t1, 
		const arma::Row<arma::uword>::fixed<3> &t2){
		
		// find common nodes
		arma::Mat<arma::uword>::fixed<3,3> M;
		for(arma::uword j=0;j<3;j++)
			M.row(j) = (t2==t1(j)).eval();

		// find common elements
		const arma::Col<arma::uword> cn1 = arma::sum(M,1);
		const arma::Row<arma::uword> cn2 = arma::sum(M,0);

		// ensure that the triangles can be merged
		assert(arma::sum(cn1)==2); assert(arma::sum(cn2)==2);

		// create quad
		arma::Row<arma::uword>::fixed<4> q;
		q(0) = arma::as_scalar(t1.cols(arma::find(cn1==false))); 
		q(1) = arma::as_scalar(t1.cols(arma::find(cn1==true,1,"first")));
		q(2) = arma::as_scalar(t2.cols(arma::find(cn2==false))); 
		q(3) = arma::as_scalar(t1.cols(arma::find(cn1==true,1,"last")));

		// return quad 
		return q;
	}

	// refinement algorithm (Catmull Clark subdivision)
	// output is pure quad mesh while input can be tri or quad
	void DistMesh2D::refine(arma::Mat<arma::uword> &qr, arma::Mat<fltp> &pr, 
		const arma::Mat<arma::uword> &qt, const arma::Mat<fltp> &p){

		// typoe of mesh
		const arma::uword num_points = p.n_rows;
		const arma::uword num_elements = qt.n_rows;
		const arma::uword num_sides = qt.n_cols;

		// list of new points one point for each side 
		// and one point in the middle of the element
		arma::Mat<fltp> pnew(num_elements*(num_sides+1),2,arma::fill::zeros);

		// list of new elements
		qr.set_size(num_elements*num_sides,4);

		// split nodes
		for(arma::uword i=0;i<num_elements;i++){
			// walk over vertices
			for(arma::uword j=0;j<num_sides;j++){
				// get indices
				const arma::uword idx1 = qt(i,j); 
				const arma::uword idx2 = qt(i,(j+1)%num_sides);

				// add new point on middle of each edge
				pnew.row(i*(num_sides+1)+j) = (p.row(idx1) + p.row(idx2))/2;
			}

			// new point on middle of element
			pnew.row((i+1)*(num_sides+1)-1) = arma::mean(p.rows(qt.row(i)),0);

			// calculate shift due to element
			const arma::uword shft = num_points + i*(num_sides+1);

			// connect new elements
			for(arma::uword j=0;j<num_sides;j++){
				// calculate last index
				const arma::uword jmo = j==0 ? num_sides-1 : j-1;

				// create quadrilateral element
				const arma::Row<arma::uword>::fixed<4> quad{qt(i,j),shft+j,shft+num_sides,shft+jmo};

				// insert element in list
				qr.row(i*num_sides+j) = quad;
			}
		}

		// add points
		pr = arma::join_vert(p,pnew);
	}

	// remove unreferenced points
	void DistMesh2D::remove_unreferenced(arma::Mat<fltp>&p, arma::Mat<arma::uword>&qt){
		// get number of nodes
		const arma::uword num_nodes = p.n_rows;

		// allocate
		arma::Col<arma::uword> node_valence(num_nodes,arma::fill::zeros);

		// calculate the valence of the nodes
		for(arma::uword i=0;i<qt.n_rows;i++)
			node_valence.rows(qt.row(i))+=1;

		// reindex array
		const arma::Col<arma::uword> reindex = arma::join_vert(
			arma::Col<arma::uword>{0}, arma::cumsum(node_valence>0));

		// shed rows
		p.shed_rows(arma::find(node_valence==0));
		qt = arma::reshape(reindex(arma::vectorise(qt)),qt.n_rows,qt.n_cols);
	}


	// remove duplicate points
	// snapping algorithm from P.O. Persson
	void DistMesh2D::remove_duplicate(){

		// check points
		if(!p_.is_finite())rat_throw_line("points contain nan or infinite");
		if(!qt_.is_finite())rat_throw_line("elements contain nan or infinite");

		// counters
		const arma::uword num_points = p_.n_rows;
		const arma::uword num_elements = qt_.n_rows;
		const arma::uword num_sides = qt_.n_cols;

		// find typical size
		const arma::Row<fltp>::fixed<2> delta = arma::max(p_,0)-arma::min(p_,0);
		const fltp dp = arma::max(delta);
		const fltp snap = dp*1024*1024*arma::Datum<fltp>::eps;

		// snap points to integers
		arma::Mat<fltp> sp = arma::round(p_/snap)*snap;

		// sort rows
		arma::Col<arma::uword> original_index = arma::regspace<arma::Col<arma::uword> >(0,num_points-1);
		for(arma::uword i=0;i<sp.n_cols;i++){
			// create sort index
			const arma::Col<arma::uword> sort_index = arma::stable_sort_index(sp.col(i));
			
			// perform sorting 
			sp = sp.rows(sort_index); original_index = original_index.rows(sort_index);
		}
		
		// find unique 
		arma::Col<arma::uword> connect = arma::regspace<arma::Col<arma::uword> >(0,num_points-1);
		arma::Col<arma::uword> duplicate(num_points,arma::fill::zeros);
		for(arma::uword i=1;i<sp.n_rows;i++){
			// check if points are identical
			if(arma::all(sp.row(i)==sp.row(i-1))){
				// connect to last non duplicate
				connect(original_index(i)) = connect(original_index(i-1));

				// mark duplicate
				duplicate(original_index(i)) = true;
			}
		}

		// shift indexes for node removal
		const arma::Col<arma::uword> idx = arma::join_vert(
			arma::Col<arma::uword>{0},
			arma::cumsum(duplicate==false));
		connect = idx(connect);

		// Remove duplicate nodes 
		p_ = p_.rows(arma::find(duplicate==false));
		qt_ = arma::reshape(connect(arma::vectorise(qt_)),num_elements,num_sides);

		{
			cmn::ShGmshFilePr gmsh = cmn::GmshFile::create("test1.gmsh");
			gmsh->write_nodes(p_.t());
			gmsh->write_elements(qt_.t());
		}

		// check
		#ifndef NDEBUG
		for(arma::uword i=0;i<qt_.n_cols;i++)assert(arma::all(qt_.col(i)!=qt_.col((i+1)%qt_.n_cols)));
		#endif
	}

	// clockwise fix
	void DistMesh2D::fix_clockwise(){
		// count
		const arma::uword num_sides = qt_.n_cols;
		const arma::uword num_elements = qt_.n_rows;

		// check if elements are clockwise
		for(arma::uword i=0;i<num_elements;i++){
			// get coordinates
			const arma::Mat<fltp> pe = p_.rows(qt_.row(i));

			// clockwise check
			fltp v = RAT_CONST(0.0);
			for(arma::uword j=0;j<num_sides;j++){
				// get indices
				const arma::uword idx1 = j, idx2 = (j+1)%num_sides;

				// calculate
				v += (pe(idx2,0) - pe(idx1,0))*(pe(idx2,1) + pe(idx1,1));
			}

			// fix element
			if(v<RAT_CONST(0.0))qt_.row(i) = arma::fliplr(qt_.row(i));
		}
	}

	// create node graph
	void DistMesh2D::create_edges(){
		// check for empty
		if(qt_.is_empty())rat_throw_line("there is no elements");

		// count
		const arma::uword num_elements = qt_.n_rows;
		const arma::uword num_sides = qt_.n_cols;

		// allocate edges
		e_.set_size(num_elements*num_sides,2);

		// copy edges
		for(arma::uword i=0;i<num_sides;i++)
			e_.rows(i*num_elements,(i+1)*num_elements-1) = 
				qt_.cols(arma::Row<arma::uword>{i,(i+1)%num_sides});

		// sort columns this means that the edge will always 
		// run from the lower index to the higher index
		e_ = arma::sort(e_,"ascend",1);

		// sort rows
		for(arma::uword i=0;i<e_.n_cols;i++){
			const arma::Col<arma::uword> sort_index = 
				arma::stable_sort_index(e_.col(i));
			e_ = e_.rows(sort_index);
		} 

		// find unique 
		arma::Col<arma::uword> is_unique(e_.n_rows,arma::fill::zeros); 
		is_unique(0) = 1; arma::uword last_unique = 0;
		for(arma::uword i=1;i<e_.n_rows;i++){
			if(arma::any(e_.row(i)!=e_.row(i-1))){
				is_unique(i)++;	last_unique = i;
			}else{
				is_unique(last_unique)++;
			}
		}

		// find surface edges
		es_ = e_.rows(arma::find(is_unique==1));

		// find internal edges
		ei_ = e_.rows(arma::find(is_unique>1));

		// select only unique rows
		e_ = e_.rows(arma::find(is_unique>0));

		// find surface node indexes
		ps_ = arma::unique(arma::vectorise(es_));

		// get only unique edges
		// e_ = cmn::Extra::unique_rows();

		// check if there are edges connected to same node
		assert(arma::all(e_.col(0)!=e_.col(1)));
	}

	// remove high valence nodes
	void DistMesh2D::high_valence_removal(const arma::uword valence){
		// check
		assert(qt_.max()<p_.n_rows);

		// check if there are edges connected to same node
		assert(arma::all(e_.col(0)!=e_.col(1)));

		// get number of nodes
		const arma::uword num_nodes = p_.n_rows;

		// allocate
		arma::Col<arma::uword> node_valence(num_nodes,arma::fill::zeros);

		// calculate the valence of the nodes
		for(arma::uword i=0;i<e_.n_rows;i++)
			node_valence.rows(e_.row(i))+=1;

		// check
		if(arma::any(node_valence==0))rat_throw_line("zero valence nodes detected.");

		// find nodes with valence of 6
		arma::Col<arma::uword> v6id = arma::find(node_valence==valence);

		// high valence node can not be on surface
		v6id.shed_rows(arma::find(cmn::Extra::ismember(v6id.t(), ps_.t())));

		// no nodes with this valence number
		if(v6id.empty())return;

		// half valence
		const arma::uword half_valence = arma::uword(std::round(fltp(valence)/2));

		// new nodes and elements
		arma::Mat<fltp> pnew(2*v6id.n_elem,2);
		arma::Mat<arma::uword> qnew(v6id.n_elem*(valence+1),qt_.n_cols);

		// create counters
		arma::uword qcnt = 0, pcnt = 0;

		// elements to be removed
		arma::Col<arma::uword> qrem(qt_.n_rows,arma::fill::zeros);

		// check which elements are already used
		arma::Col<arma::uword> is_used(qt_.n_rows,arma::fill::zeros);

		// walk over high valence nodes
		for(arma::uword i=0;i<v6id.n_elem;i++){
			// get index of high valence node from list
			const arma::uword hvn = v6id(i);

			// get neighbouring elements
			const arma::Col<arma::uword> nid = arma::find(arma::any(qt_==hvn,1));
			assert(nid.n_elem==valence);

			// check if the elements are already used
			if(arma::any(is_used(nid)))continue;

			// mark nodes used
			is_used(nid).fill(1);

			// get local mesh to be replaced
			arma::Mat<arma::uword> qt_local = qt_.rows(nid);
			assert(qt_local.n_rows==valence);

			// rotate nodes
			for(arma::uword j=0;j<valence;j++){
				// find high valence node in this row
				const arma::Col<arma::uword> list = arma::find(qt_local.row(j)==hvn);
				assert(!list.empty());

				// get position of the valence node
				const arma::uword pos = arma::as_scalar(list);

				// shift valence node to first position
				qt_local.row(j) = arma::shift(qt_local.row(j),-arma::sword(pos));

			}
			assert(arma::all(qt_local.col(0)==qt_local(0,0)));

			// find overlapping nodes and add them to the list
			const arma::Mat<arma::uword> list = qt_local.cols(arma::Col<arma::uword>{1,3});
			assert(arma::unique(list.col(0)).eval().n_elem==valence);
			assert(arma::unique(list.col(1)).eval().n_elem==valence);

			// link list
			arma::Col<arma::uword> linkedlist(valence); linkedlist(0) = 0;
			for(arma::uword j=1;j<valence;j++)
				linkedlist(j) = arma::as_scalar(arma::find(list(linkedlist(j-1),1)==list.col(0)));

			// find clockwise mesh
			arma::Row<fltp> pv = p_.row(hvn);
			arma::Row<fltp> p1 = p_.row(qt_local(linkedlist(1),2));
			arma::Row<fltp> p2 = p_.row(qt_local(linkedlist(1+half_valence),2));

			// compute two new vertices
			arma::Row<fltp> pv1 = (2*pv + p1)/3;
			arma::Row<fltp> pv2 = (2*pv + p2)/3;

			// add to new node table
			pnew.row(pcnt) = pv1;
			pnew.row(pcnt+1) = pv2;

			// index of these new nodes
			const arma::uword vnidx1 = num_nodes+pcnt;
			const arma::uword vnidx2 = num_nodes+pcnt+1;

			// modify EToV table and insert new element
			for(arma::uword j=0;j<half_valence;j++)qt_local(linkedlist(j),0) = vnidx1;
			for(arma::uword j=half_valence;j<valence;j++)qt_local(linkedlist(j),0) = vnidx2;

			// add new element (always split by 1st element)
			const arma::Row<arma::uword> qadd{qt_local(linkedlist(0),1),vnidx1,qt_local(linkedlist(half_valence-1),3),vnidx2};

			// add face
			qt_local = arma::join_vert(qt_local, qadd);

			// add to new elements
			qnew.rows(qcnt,qcnt+valence+1-1) = qt_local;

			// mark old elements for removal
			qrem.rows(nid).fill(1);

			// increment counters
			qcnt += valence+1; pcnt += 2;
		}

		// add new points and new elements to old points and elements
		qt_ = arma::join_vert(qt_,qnew.rows(0,qcnt-1));
		p_ = arma::join_vert(p_,pnew.rows(0,pcnt-1));

		// check
		assert(qt_.max()<p_.n_rows);

		// create re-indexing
		arma::Col<arma::uword> pkeep(p_.n_rows,arma::fill::ones); pkeep(v6id).fill(0);
		arma::Col<arma::uword> reidx(p_.n_rows); arma::uword cnt=0;
		for(arma::uword i=0;i<p_.n_rows;i++){
			reidx(i) = cnt; cnt+=pkeep(i);
		}

		// drop points
		p_.shed_rows(v6id);

		// drop elements
		qt_.shed_rows(arma::find(qrem));

		// update indices
		qt_ = arma::reshape(reidx(arma::vectorise(qt_)),qt_.n_rows,qt_.n_cols);

		// check
		assert(qt_.max()<p_.n_rows);

		// re-create edge list
		create_edges();
	}

	// recombination from triangles to quad
	void DistMesh2D::tri2quad(){
		// check tri mesh
		if(qt_.n_cols!=3)rat_throw_line("elements are not triangular");

		// check if there are edges connected to same node
		assert(arma::all(e_.col(0)!=e_.col(1)));

		// check mesh
		assert(arma::all(cmn::Triangle::calc_area(arma::join_vert(p_.t(),arma::Row<fltp>(p_.n_rows,arma::fill::zeros)),qt_.t())>0));

		// get element graph
		arma::Col<arma::uword> row,col;
		element_graph(row,col,qt_);

		// allocate all possible quads
		arma::Mat<arma::uword> qnew(row.n_elem,4);

		// walk over combinations and calculate element quality
		for(arma::uword i=0;i<row.n_elem;i++){	
			// find common nodes
			const arma::Row<arma::uword>::fixed<3> t1 = qt_.row(row(i));
			const arma::Row<arma::uword>::fixed<3> t2 = qt_.row(col(i));

			// merge
			qnew.row(i) = merge2triangles(t1,t2);
		}

		// get quality for proposed elements
		arma::Col<fltp> Q;
		if(alternate_recombination_){
			Q = element2quality_alt(qnew,p_);
		}else{
			Q = element2quality(qnew,p_);
		}

		// greedy greedy algorithm ...
		// selects best quads for itself and 
		// leaves bad triangles for catmull clark subdivision
		// sort idexes. This could possibly be replaced with a perfect 
		// matching algorithm.
		arma::Col<arma::uword> srt = arma::sort_index(Q,"descend");

		// create original index
		arma::Col<arma::uword> original_index = arma::regspace<arma::Col<arma::uword> >(0,Q.n_rows-1);
		original_index = original_index.rows(srt);

		// keep track of triangles used
		arma::Col<arma::uword> triangle_used(qt_.n_rows,arma::fill::zeros);
		arma::Col<arma::uword> quad_used(qnew.n_rows,arma::fill::zeros);

		// select
		for(arma::uword i=0;i<row.n_elem;i++){
			// check quality
			if(Q(original_index(i))>recombination_treshold_){
				// get indexes
				arma::uword qidx = original_index(i);
				arma::uword tidx1 = row(qidx);
				arma::uword tidx2 = col(qidx);

				// check if both triangles still available
				if(triangle_used(tidx1)==false && triangle_used(tidx2)==false){
					// mark quad for creation
					quad_used(qidx) = true;
					triangle_used(tidx1) = true; 
					triangle_used(tidx2) = true;
				}
			}
		}

		// get new quads from list
		qnew = qnew.rows(arma::find(quad_used));

		// get unused triangles from list
		arma::Mat<arma::uword> tremain = qt_.rows(arma::find(triangle_used==false));

		// perform catmul clark on triangles and split rectangles
		arma::Mat<arma::uword> q1, q2; arma::Mat<fltp> p1, p2;
		refine(q1,p1,qnew,p_); refine(q2,p2,tremain,p_);

		// merge meshes
		qt_ = arma::join_vert(q1,q2+p1.n_rows);
		p_ = arma::join_vert(p1,p2);

		// remove duplicate nodes
		remove_duplicate();

		// fix clockwise
		fix_clockwise();

		// re-create edge list
		create_edges();

		// find snapping size
		fltp snap = arma::as_scalar(
			arma::max(arma::max(bbox_,0)-arma::min(bbox_,0),1))*
			1024*arma::Datum<fltp>::eps;

		// find points on surface that are not fixed
		arma::Col<arma::uword> is_unique(ps_.n_rows,arma::fill::ones);
		for(arma::uword i=0;i<pfix_.n_rows;i++){
			for(arma::uword j=0;j<ps_.n_rows;j++){
				is_unique(j) &= !arma::all(
					arma::round(pfix_.row(i)/snap)*snap == 
					arma::round(p_.row(ps_(j))/snap)*snap);
			}
		}
		arma::Col<arma::uword> psnf = ps_(arma::find(is_unique));

		// calculate d for surface points
		arma::Mat<fltp> ps = p_.rows(psnf);
		arma::Col<fltp> ds = fd_->calc_distance(ps);

		// calculate numerical gradient of the distance function
		arma::Mat<fltp> pdx = p_.rows(psnf); pdx.col(0)+=deps_;
		arma::Mat<fltp> pdy = p_.rows(psnf); pdy.col(1)+=deps_;
		arma::Col<fltp> dgradx = (fd_->calc_distance(pdx)-ds)/deps_;
		arma::Col<fltp> dgrady = (fd_->calc_distance(pdy)-ds)/deps_;
		arma::Col<fltp> dgrad2 = dgradx%dgradx + dgrady%dgrady;

		// Project points back to boundary
		p_.rows(psnf) -= arma::join_horiz(ds%dgradx/dgrad2, ds%dgrady/dgrad2);

		// fix clockwise
		fix_clockwise();

		// remove high valence nodes
		high_valence_removal(6);
		high_valence_removal(6);

		// smooth mesh
		smoothen();
	}


	// get node valence
	arma::Col<arma::uword> DistMesh2D::get_node_valence()const{
		// allocate
		arma::Col<arma::uword> node_valence(p_.n_rows,arma::fill::zeros);

		// calculate the valence of the nodes
		for(arma::uword i=0;i<e_.n_rows;i++)
			node_valence.rows(e_.row(i))+=1;

		return node_valence;
	}

	// smooth triangle mesh
	void DistMesh2D::smoothen(){
		// settings
		const arma::uword num_iter = 10;
		const fltp dmp = RAT_CONST(0.7);

		// create neighbour list
		arma::Mat<arma::uword> nb = arma::join_vert(e_,arma::fliplr(e_));
		arma::Col<arma::uword> sort_index = arma::sort_index(nb.col(0));
		nb = nb.rows(sort_index);

		// find sections
		arma::Mat<arma::uword> sects = cmn::Extra::find_sections(nb.col(0).t());

		// iterations
		for(arma::uword j=0;j<num_iter;j++){
			// copy points list
			arma::Mat<fltp> pnew = p_;

			// walk over sections
			for(arma::uword i=0;i<sects.n_cols;i++){
				// node index
				arma::uword myidx = nb(sects(0,i),0);

				// find neighbours of this node
				arma::Col<arma::uword> mynb = nb.submat(sects(0,i),1,sects(1,i),1);

				// calculate new position
				pnew.row(myidx) = arma::mean(p_.rows(mynb),0);
			}

			// do not move surface elements
			pnew.rows(ps_) = p_.rows(ps_);

			// apply damping and move points
			p_ += dmp*(pnew-p_); 
		}
	}


	// set element size
	void DistMesh2D::set_h0(const fltp h0){
		h0_ = h0; 
		deps_ = std::sqrt(RAT_CONST(1e-7))*h0;
		geps_ = RAT_CONST(0.001)*h0;
	}

	// set seed for rng
	void DistMesh2D::set_rng_seed(const unsigned int rng_seed){
		rng_seed_ = rng_seed;
	}

	// set strength for rng
	void DistMesh2D::set_rng_strength(const fltp rng_strength){
		rng_strength_ = rng_strength;
	}

	// set quad mesh
	void DistMesh2D::set_quad(const bool quad){
		quad_ = quad;
	}

	// set distance function
	void DistMesh2D::set_distfun(const ShDistFunPr &fd){
		fd_ = fd;
	}

	// set scaling function
	void DistMesh2D::set_scalefun(const ShDistFunPr &fh){
		fh_ = fh;
	}

	// set maximum number of iterations
	void DistMesh2D::set_maxiter(const arma::uword maxiter){
		maxiter_ = maxiter;
	}

	// set fixed points
	void DistMesh2D::set_fixed(const arma::Mat<fltp> &pfix_ext){
		if(pfix_ext.n_cols!=2)rat_throw_line("fixed points must be supplied in 2 column matrix");
		pfix_ext_ = pfix_ext;
	}

	// set number of node limit
	void DistMesh2D::set_num_node_limit(const arma::uword num_node_limit){
		num_node_limit_ = num_node_limit;
	}

	// set bounding box
	// void DistMesh2D::set_bounding(const arma::Mat<fltp> &bbox){
	// 	bbox_ = bbox;
	// }

	// get elements of mesh
	const arma::Mat<arma::uword>& DistMesh2D::get_elements() const{
		return qt_;
	}

	// get points of mesh
	const arma::Mat<fltp>& DistMesh2D::get_nodes() const{
		return p_;
	}

	// get surface edges
	const arma::Mat<arma::uword>& DistMesh2D::get_surface_edges() const{
		return es_;
	}

	// get surface edges
	const arma::Mat<arma::uword>& DistMesh2D::get_internal_edges() const{
		return ei_;
	}

	// get surface node indexes
	const arma::Col<arma::uword>& DistMesh2D::get_surface_node_idx() const{
		return ps_;
	}

	const ShDistFunPr& DistMesh2D::get_distfun()const{
		return fd_;
	}

	const ShDistFunPr& DistMesh2D::get_scalefun()const{
		return fh_;
	}

	// distmesh setup function
	// translated from distmesh from Per Olof Persson
	// http://persson.berkeley.edu/distmesh/
	void DistMesh2D::setup(){
		// check input
		if(h0_<=0)rat_throw_line("element size must be larger than zero");
		if(pfix_ext_.n_cols!=2)rat_throw_line("fixed points must be stored in two columns x and y");

		// create a random number generator
		// we use this instead of armadillo generator to ensure each thread
		// has a reproduceable mesh without race conditions 
		std::mt19937 generator(rng_seed_); // Standard mersenne_twister_engine seeded with seed
		std::uniform_real_distribution<fltp> distribution(0.0, 1.0);

		// get bounding box
		bbox_ = fd_->get_bounding();

		// check bounding
		if(arma::any(bbox_.row(0)>=bbox_.row(1)))rat_throw_line("negatively or zero spaced bounding box");

		// extend bounding
		bbox_.row(0) -= RAT_CONST(0.02)*(bbox_.row(1) - bbox_.row(0));
		bbox_.row(1) += RAT_CONST(0.02)*(bbox_.row(1) - bbox_.row(0));

		// check if bounding box has any volume
		if(arma::prod(bbox_.row(1)-bbox_.row(0))<=0)rat_throw_line("bounding box inconsistent");

		// get fixed points
		const fltp abstol = h0_/20;
		pfix_ = arma::join_vert(fd_->get_fixed(abstol),pfix_ext_);

		// remove duplicate pfix points
		arma::Col<arma::uword> mark_duplicate_pfix(pfix_.n_rows,arma::fill::zeros);
		for(arma::uword i=0;i<pfix_.n_rows;i++){
			for(arma::uword j=i+1;j<pfix_.n_rows;j++){
				if(arma::all(arma::abs(pfix_.row(i) - pfix_.row(j))<h0_/20))
					mark_duplicate_pfix(j) = 1;
			}
		}
		pfix_.shed_rows(arma::find(mark_duplicate_pfix));

		// check fixed points
		assert(pfix_.is_finite());

		// check points
		if(pfix_.n_cols!=2)rat_throw_line("fixed points must be supplied in 2 column matrix");

		// check which fixed points are on the boundary and discard others
		pfix_ = pfix_.rows(arma::find(arma::abs(fd_->calc_distance(pfix_))<abstol));

		// check number of elements
		if(((bbox_(1,0) - bbox_(0,0))/h0_)*((bbox_(1,1) - bbox_(0,1))/(h0_*std::sqrt(3)/2))>num_node_limit_)
			rat_throw_line("node limit exceeded, if needed it can be increased or increase the element size");

		// 1. create initial distribution in bounding box (equilateral triangles)
		// generate grid of coordinates
		arma::Row<fltp> xgv = arma::regspace<arma::Row<fltp> >(bbox_(0,0),h0_,bbox_(1,0));
		arma::Col<fltp> ygv = arma::regspace<arma::Col<fltp> >(bbox_(0,1),h0_*std::sqrt(3)/2,bbox_(1,1));
		arma::Mat<fltp> x,y;
		meshgrid(x,y,xgv,ygv);

		// shift even rows
		x.rows(arma::regspace<arma::Col<arma::uword> >(2,2,x.n_rows-1)) += h0_/2;

		// create list of node coordinates
		p_ = arma::join_horiz(arma::vectorise(x),arma::vectorise(y));

		// add some random noise on the points
		for(arma::uword i=0;i<p_.n_rows;i++){
			p_(i,0) += (2*distribution(generator) - 1.0)*rng_strength_*h0_;
			p_(i,1) += (2*distribution(generator) - 1.0)*rng_strength_*h0_;
		}

		// check node coordinates
		assert(p_.is_finite());

		// 2. remove points outside the region, apply the rejection method
		// evaluate distance function
		arma::Col<fltp> d = fd_->calc_distance(p_);

		// keep only d<0 points
		p_ = p_.rows(arma::find(d<geps_));

		// check if any point left
		if(p_.is_empty())rat_throw_line("shape function inconsistent");

		// probability to keep point
		arma::Col<fltp> r0 = RAT_CONST(1.0)/arma::pow(arma::clamp(fh_->calc_distance(p_),h0_,arma::Datum<fltp>::inf),2);

		// create random number for rejection method
		arma::Col<fltp> randu(p_.n_rows);
		for(arma::uword i=0;i<p_.n_rows;i++)randu(i) = distribution(generator);

		// rejection method
		p_ = p_.rows(arma::find(randu<r0/arma::max(r0)));

		// find snapping size
		const fltp snap = arma::as_scalar(arma::max(arma::max(bbox_,0)-arma::min(bbox_,0),1))*1024*arma::Datum<fltp>::eps;

		// Remove duplicated nodes
		arma::Col<arma::uword> is_unique(p_.n_rows,arma::fill::ones);
		for(arma::uword i=0;i<pfix_.n_rows;i++){
			for(arma::uword j=0;j<p_.n_rows;j++){
				is_unique(j) &= !arma::all(
					arma::round(pfix_.row(i)/snap)*snap == 
					arma::round(p_.row(j)/snap)*snap);
			}
		}
		p_ = p_.rows(arma::find(is_unique));
		arma::uword num_fix = pfix_.n_rows;

		// check number of points
		if(p_.is_empty())rat_throw_line("all points were removed");

		// prepend fix points
		p_ = arma::join_vert(pfix_,p_);

		// check node coordinates
		assert(p_.is_finite());

		// number of points N
		arma::uword N = p_.n_rows;

		// setup for first iteration
		arma::Mat<fltp> pold(N,2); pold.fill(arma::Datum<fltp>::inf);

		// iterate
		for(arma::uword k=0;k<maxiter_;k++){
			// check pold
			assert(!pold.is_empty());

			// check node coordinates
			assert(p_.is_finite());

			// count++;
			// 3. retriangulation by the Delaunay algorithm
			// Any large movement?
			//if(arma::as_scalar(arma::any(arma::sqrt(arma::sum(arma::pow(p_-pold,2),1))>ttol))){
			fltp delta = arma::as_scalar(arma::max(arma::sqrt(arma::sum(arma::pow(p_-pold,2),1))/h0_));
			if(delta>ttol_){
				// save current positions
				pold = p_;                                          

				// create list of triangles
				qt_ = triangulation(p_); 

				// Compute centroids
				const arma::Mat<fltp> pmid = (p_.rows(qt_.col(0)) + p_.rows(qt_.col(1)) + p_.rows(qt_.col(2)))/3;

				// Keep interior triangles only
				qt_ = qt_.rows(arma::find(fd_->calc_distance(pmid)<(-geps_)));
				// arma::uword num_tri = qt_.n_rows;
				if(qt_.n_rows==0)rat_throw_line("no triangles left inside mesh");

				// create edge list
				create_edges();

				// // 4. Describe each bar by a unique pair of nodes
				// // Interior bars duplicated
				// bars.set_size(3*num_tri,2);
				// bars.rows(0*num_tri,1*num_tri-1) = qt_.cols(arma::Row<arma::uword>{0,1});
				// bars.rows(1*num_tri,2*num_tri-1) = qt_.cols(arma::Row<arma::uword>{0,2});
				// bars.rows(2*num_tri,3*num_tri-1) = qt_.cols(arma::Row<arma::uword>{1,2});

				// // get only unique edges
				// bars = cmn::Extra::unique_rows(arma::sort(bars,"ascend",1));
			}

			// 6. Move mesh points based on bar lengths L and forces F
			// List of bar vectors
			const arma::Mat<fltp> evec = p_.rows(e_.col(0))-p_.rows(e_.col(1));

			// L = Bar lengths
			arma::Col<fltp> L = arma::sqrt(arma::sum(evec%evec,1));

			// Scaling function
			arma::Mat<fltp> pbar = (p_.rows(e_.col(0))+p_.rows(e_.col(1)))/2;
			arma::Col<fltp> he = arma::clamp(fh_->calc_distance(pbar),h0_,arma::Datum<fltp>::inf);

			// L0 = Desired lengths
			arma::Col<fltp> L0 = he*Fscale_*std::sqrt(arma::sum(L%L)/arma::sum(he%he));     

			// Density control - remove points that are too close
			if((k+1)%densityctrlfreq_==0 && arma::any(L0>2*L)){
				// std::cout<<"controlling density"<<std::endl;
				// get points to remove
				arma::Col<arma::uword> idx_remove = arma::vectorise(e_.rows(arma::find(L0>2*L)));

				// do not touch fixed
				idx_remove = idx_remove.rows(arma::find(idx_remove>=num_fix));

				// remove rows from p
				arma::Col<arma::uword> idx_keep(N,arma::fill::ones);
				idx_keep.rows(idx_remove).fill(0);
				p_ = p_.rows(arma::find(idx_keep));

				// update counter
				N = p_.n_rows;	

				// ensure that the mesh is re-setup	
				pold.set_size(N,2);	pold.fill(arma::Datum<fltp>::inf);

				// check if there are points remaining
				if(p_.is_empty())rat_throw_line("all points were removed during density control");

				// restart
				continue;
			}

			// check node coordinates
			assert(p_.is_finite());

			// calculate bar forces (scalars)
			arma::Col<fltp> F = arma::clamp(L0-L,0,arma::Datum<fltp>::inf); 
			
			// Bar forces (x,y components)
			// arma::Mat<fltp> Fvec = (F/L)%evec.each_col();
			arma::Mat<fltp> Fvec = F/L*arma::Row<fltp>{1,1}%evec;

			// Force on points
			arma::Mat<fltp> Ftot(N,2,arma::fill::zeros);
			for(arma::uword j=0;j<e_.n_rows;j++){
				Ftot.row(e_(j,0)) += Fvec.row(j);
				Ftot.row(e_(j,1)) -= Fvec.row(j);
			}
			// Ftot=full(sparse(bars(:,[1,1,2,2]),ones(size(F))*[1,2,1,2],[Fvec,-Fvec],N,2));

			// Force = 0 at fixed points
			if(num_fix>0)Ftot.rows(0,num_fix-1).fill(0); 
			
			// Update node positions
			p_ += deltat_*Ftot;

			// 7. Bring outside points back to the boundary
			// Find points outside (d>0)
			d=fd_->calc_distance(p_); assert(!d.is_empty());
			arma::Col<arma::uword> idx_outside = arma::find(d>0);

			// check if any internal points left
			// assert(arma::any(d<(-geps_)));

			// calculate numerical gradient of the distance function
			arma::Mat<fltp> pdx = p_.rows(idx_outside); pdx.col(0)+=deps_;
			arma::Mat<fltp> pdy = p_.rows(idx_outside); pdy.col(1)+=deps_;
			arma::Col<fltp> dgradx = (fd_->calc_distance(pdx)-d.rows(idx_outside))/deps_;
			arma::Col<fltp> dgrady = (fd_->calc_distance(pdy)-d.rows(idx_outside))/deps_;
			// arma::Col<fltp> dgrad2 = dgradx%dgradx + dgrady%dgrady;
			arma::Col<fltp> dgrad2 = arma::hypot(dgradx,dgrady); // this is not the same as the original shown in the previous commented line

			// check gradients
			assert(dgradx.is_finite()); assert(dgrady.is_finite());

			// gradient can not be zero
			dgrad2(arma::find(dgrad2==0)).fill(1.0);

			// Project points back to boundary
			p_.rows(idx_outside) -= arma::join_horiz(d(idx_outside)%dgradx/dgrad2, d(idx_outside)%dgrady/dgrad2);

			// check node coordinates
			assert(p_.is_finite());

			// drop elements that now have zero area
			qt_.shed_rows(arma::find(cmn::Triangle::calc_area(arma::join_vert(p_.t(),arma::Row<fltp>(p_.n_rows,arma::fill::zeros)),qt_.t())<=1e-9));

			// 8. Termination criterion: All interior nodes move less than dptol (scaled)
			const arma::Col<arma::uword> idx = arma::find(d<(-geps_));
			if(!idx.is_empty()){
				fltp conv = arma::as_scalar(arma::max(arma::sqrt(arma::sum(
					deltat_*arma::square(Ftot.rows(idx)),1))/h0_));

				// std::cout<<conv<<" "<<dptol<<std::endl;
				if(conv<dptol_ && k!=0)break;
			}

			// if max(sqrt(sum(deltat*Ftot(d<-geps,:).^2,2))/h0)<dptol, break; end
			//max(sqrt(sum(deltat*Ftot(d<-geps,:).^2,2))/h0)<dptol, break; end
		}

		// remove unreferenced points
		remove_unreferenced(p_,qt_);

		// remove duplicate nodes (if any)
		remove_duplicate();

		// fix clockwise (if needed)
		fix_clockwise();

		// Quad mesh recombination
		if(quad_)tri2quad();
	}

	// gmsh export
	void DistMesh2D::export_gmsh(const cmn::ShGmshFilePr& gmsh) const{
		// write points in xy plane with z=0
		gmsh->write_nodes(p_.t());

		// write elements
		gmsh->write_elements(qt_.t());
	}

	// gmsh export
	void DistMesh2D::export_gmsh_perimeter(const cmn::ShGmshFilePr& gmsh) const{
		// create perimeter
		const ShPerimeterPr perimeter = fd_->create_perimeter(fd_->calc_bounding_size()/64);

		// write points in xy plane with z=0
		gmsh->write_nodes(perimeter->get_nodes());

		// write elements
		gmsh->write_elements(perimeter->get_elements());
	}


	// function for calculating area
	fltp DistMesh2D::calc_area() const{
		// get counters
		const arma::uword num_elements = qt_.n_rows;
		const arma::uword num_edges = qt_.n_cols;

		// allocate areas
		arma::Col<fltp> A(num_elements,arma::fill::zeros);

		// calculate areas
		for(arma::uword i=0;i<num_elements;i++){
			// get points
			const arma::Mat<fltp> pt = p_.rows(qt_.row(i));

			// walk over triangles or quadrangles
			for(arma::uword j=0;j<num_edges-2;j++){
				// get two edges
				const arma::Row<fltp>::fixed<2> a = pt.row((2*j+1)%num_edges) - pt.row((2*j)%num_edges);
				const arma::Row<fltp>::fixed<2> b = pt.row((2*j+2)%num_edges) - pt.row((2*j+1)%num_edges);

				// calculate area using cross product
				A(i) += std::abs(a(0)*b(1) - a(1)*b(0))/2;
			}
		}

		// sum areas and return
		return arma::sum(A);
	}

	// algorithm for making 3D meshes
	void DistMesh2D::extrude(const fltp ell, const arma::uword num_z){
		// create z coords
		const arma::Col<fltp> z = arma::linspace<arma::Col<fltp> >(-ell/2,ell/2,num_z+1);
		
		const arma::uword num_elements_2d = qt_.n_rows;
		const arma::uword num_nodes_2d = p_.n_rows;

		// copy nodes into layers
		arma::Mat<fltp> p3d((num_z+1)*num_nodes_2d, 3);
		for(arma::uword i=0;i<num_z+1;i++){
			p3d.rows(i*num_nodes_2d,(i+1)*num_nodes_2d-1) = 
				arma::join_horiz(p_,arma::Col<fltp>(num_nodes_2d,arma::fill::value(z(i))));
		}

		// copy elements
		arma::Mat<arma::uword> h(num_z*num_elements_2d, 2*qt_.n_cols);
		for(arma::uword i=0;i<num_z;i++){
			h.rows(i*num_elements_2d,(i+1)*num_elements_2d-1) = 
				arma::join_horiz(qt_+ i*num_nodes_2d, qt_+ (i+1)*num_nodes_2d);
		}

		// set to self
		qt_ = h; p_ = p3d;
	}

	// type string for serialization
	std::string DistMesh2D::get_type(){
		return "rat::dm::dfdistmesh2d";
	}

	// method for serialization into json
	void DistMesh2D::serialize(Json::Value &js, cmn::SList &list) const{
		// parent
		cmn::Node::serialize(js,list);

		// type
		js["type"] = get_type();

		// quad mesher enabled
		js["quad"] = quad_;

		// element size
		js["h0"] = h0_;
		
		// externally supplied fixed point list
		for(arma::uword i=0;i<pfix_ext_.n_rows;i++){
			js["pfix_x"].append(pfix_ext_(i,0));
			js["pfix_y"].append(pfix_ext_(i,1));
		}

		// distance functions
		js["geometryfun"] = cmn::Node::serialize_node(fd_,list);
		js["scalingfun"] = cmn::Node::serialize_node(fh_,list);
		js["num_node_limit"] = static_cast<int>(num_node_limit_);

		// random number generator
		js["rng_seed"] = rng_seed_;
		js["rng_strength"] = rng_strength_;
	}

	// method for deserialisation from json
	void DistMesh2D::deserialize(
		const Json::Value &js, cmn::DSList &list, 
		const cmn::NodeFactoryMap &factory_list, 
		const boost::filesystem::path &pth){
		
		// parent
		cmn::Node::deserialize(js,list,factory_list,pth);

		// quad mesher enabled 
		set_quad(js["quad"].asBool());

		// element size
		set_h0(js["h0"].ASFLTP());

		// externally supplied fixed point list
		pfix_ext_.set_size(js["pfix_x"].size(),2);
		arma::uword idx = 0; auto it1 = js["pfix_x"].begin(); auto it2 = js["pfix_y"].begin();
		for(;it1!=js["pfix_x"].end() && it2!=js["pfix_y"].end();it1++,it2++,idx++){
			pfix_ext_(idx,0) = (*it1).ASFLTP(); 
			pfix_ext_(idx,1) = (*it2).ASFLTP();
		}
		
		// distance functions
		fd_ = cmn::Node::deserialize_node<DistFun>(js["geometryfun"], list, factory_list, pth);
		fh_ = cmn::Node::deserialize_node<DistFun>(js["scalingfun"], list, factory_list, pth);

		// node limit
		set_num_node_limit(js["num_node_limit"].asUInt64());

		// random number generator
		rng_seed_ = js["rng_seed"].asUInt();
		rng_strength_ = js["rng_strength"].asUInt();
	}

	// // custom function for finding sections in array
	// // for example find_sections([1,1,1,2,2,2,3,3],"first") should give [0,3,6]
	// // for example find_sections([1,1,1,2,2,2,3,3],"last") should give [2,5,7]
	// arma::Mat<arma::uword> DistMesh2D::find_sections(const arma::Row<arma::uword> &M){
	// 	// make sure it is a row vector
	// 	assert(M.n_rows==1);

	// 	// allocate output
	// 	arma::Mat<arma::uword> indices;

	// 	// check if matrix is empty
	// 	if(!M.is_empty()){
	// 		// find indices at changing parts
	// 		arma::Row<arma::uword> idx = arma::find(M.tail_cols(M.n_cols-1)!=M.head_cols(M.n_cols-1)).t(); 
		
	// 		// make first and last indexes
	// 		arma::Row<arma::uword> idx1 = arma::join_horiz(arma::Mat<arma::uword>(1,1,arma::fill::zeros),idx+1);
	// 		arma::Row<arma::uword> idx2 = arma::join_horiz(idx,arma::Mat<arma::uword>(1,1,arma::fill::ones)*(M.n_cols-1));
		
	// 		// return
	// 		indices = arma::join_vert(idx1,idx2);
	// 	}

	// 	// if matrix is empty return empty indices
	// 	else{
	// 		indices = arma::Mat<arma::uword>(2,0);
	// 	}

	// 	// return output matrix with indexes
	// 	return indices;
	// }

}}
